%This program computes the elasticity at every .5 hour for supply and
%demand for all TORQ stocks and records the NYSE bid-ask spread quote.
%Buyside  -> OSIDE = 668589 or 668945
%Sellsise -> OSIDE = 836976 or 837643 or 838384 or 838388
clear all
load SOD
load CQ

%Set times
TIMES=[10:.5:16];

%Get list of dates from CQ
D1=sortrows(CQ(find(CQ(:,1)==706888),2));
D2=[1;D1(1:end-1)];
DATES=D1(D1~=D2);
clear D1 D2

%Get list of symbols from CQ
S1=sortrows(CQ(find(CQ(:,2)==901101),1));
S2=[1;S1(1:end-1)];
SYMBOLS=S1(S1~=S2);
clear S1 S2

%Placeholder for elasticity estimates
E=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
CON=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
TSTATCON=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
TSTATE=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
R2=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
N=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
BIDASK=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);
DEPTH=NaN*ones(length(SYMBOLS),length(DATES),length(TIMES),2);

for s=1:length(SYMBOLS)
    SYM=SYMBOLS(s);%Symbol FDX
    clear SYMIND SYMINDQT
    SYMIND=find(SOD(:,38)==SYM);
    SYMINDQT=find(CQ(:,1)==SYM);
    for d=1:length(DATES)
        DT =DATES(d); % Date
        clc
        disp(['SYM = ' num2str(s) '  DATE = ' num2str(d) ])  %TIME = ' num2str(t) '  BUY/SELL = ' num2str(b) ])
        for t=1:length(TIMES)
            %for t=11;
            TM =TIMES(t)*3600; %Time Noon
            for b=1:2
                BS =b; % Buy/Sell side indicator (Buy = 1, Sell = 2)
                clear ALLIND QTIND QT QTSIZ
                %Quote
                QTIND=find((CQ(SYMINDQT,1)==SYM)&(CQ(SYMINDQT,2)==DT)&(CQ(SYMINDQT,3)<=TM)&(CQ(SYMINDQT,10)==78));
                if isempty(QTIND)==0
                    QT=[CQ(max(SYMINDQT(QTIND)),[4 5])];
                    QTSIZ=[CQ(max(SYMINDQT(QTIND)),[6 7])];
                    if BS==1
                        ALLIND=...%All records with time and date before event time index
                        find((SOD(SYMIND,38)==SYM)...%Symbol
                        &((SOD(SYMIND,26)==49)|(SOD(SYMIND,26)==51)|(SOD(SYMIND,26)==57))...%Limit, Stop Limit, or Limit or Better 
                        &(((SOD(SYMIND,27)==DT)&(SOD(SYMIND,33)<=TM))|(SOD(SYMIND,27)<DT))...%Order date and time
                        &((SOD(SYMIND,32)==668589)|(SOD(SYMIND,32)==668945))...%OSIDE Buy
                        );
                    else
                        ALLIND=...%All records with time and date before event time index
                        find((SOD(SYMIND,38)==SYM)...%Symbol
                        &((SOD(SYMIND,26)==49)|(SOD(SYMIND,26)==51)|(SOD(SYMIND,26)==57))...%Limit, Stop Limit, or Limit or Better 
                        &(((SOD(SYMIND,27)==DT)&(SOD(SYMIND,33)<=TM))|(SOD(SYMIND,27)<DT))...%Order date and time
                        &((SOD(SYMIND,32)==836976)|(SOD(SYMIND,32)==837643)|(SOD(SYMIND,32)==838384)|(SOD(SYMIND,32)==838388))...%OSIDE Selling
                        );
                    end
                    OALL=SOD(SYMIND(ALLIND),:);
                    EIND=zeros(length(ALLIND),1);%Execution index
                    CIND=zeros(length(ALLIND),1);%Cancellation index
                    DIND=zeros(length(ALLIND),1);%Cancellation index
                    SPIND=zeros(length(ALLIND),1);%Cancellation index
                    PIND=zeros(length(ALLIND),1);%Cancellation index
                    for i=1:length(ALLIND);
                        if OALL(i,27)>=901101 %current records
                                %Executions
                            if ...%if 1 record in 1 part and executed before date time in toto
                               (OALL(i,22)==1)&&(((OALL(i,40)==DT)&&(OALL(i,36)<=TM))||(OALL(i,40)<DT))&&(OALL(i,35)>0)&&(OALL(i,30)>48)  
                                EIND(i,1)=1;
                            elseif ...%if multiple records in group and executed before date time
                               (OALL(i,22)>1)&&(((OALL(i,40)==DT)&&(OALL(i,36)<=TM))||((OALL(i,40)<DT)))&&(OALL(i,35)>0)  
                                EIND(i,1)=1;    
                            elseif ...%Original order filled later
                               (OALL(i,30)>48)&&(OALL(i,35)==0)&&(OALL(i,16)==0); 
                               %then find if there exist matches with executions or cancellations before event time
                               EIND(i,1)=isempty(find(((OALL(:,30)==48)&(OALL(:,23)==OALL(i,23))&(OALL(:,31)==OALL(i,31))&(OALL(:,38)==OALL(i,38))&(OALL(:,27)==OALL(i,27))&(OALL(:,26)==OALL(i,26))&(OALL(:,29)==OALL(i,29))&(OALL(:,32)==OALL(i,32))&(OALL(:,41)==OALL(i,41))&(((OALL(:,40)==DT)&(OALL(:,36)<=TM))|((OALL(:,40)<DT)))&(OALL(:,35)>0)), 1))==0;
                               %or if a day order goes unfilled by the end of the day
                               CIND(i,1)=(isempty(find(((OALL(:,30)==48)&(OALL(:,23)==OALL(i,23))&(OALL(:,31)==OALL(i,31))&(OALL(:,38)==OALL(i,38))&(OALL(:,27)==OALL(i,27))&(OALL(:,26)==OALL(i,26))&(OALL(:,29)==OALL(i,29))&(OALL(:,32)==OALL(i,32))&(OALL(:,41)==OALL(i,41))&(((OALL(:,40)==DT)&(OALL(:,3 )<=TM))|((OALL(:,40)<DT)))&(OALL(:,16)>0)), 1))==0)|(((OALL(i,41)==49)|(OALL(i,41)==51)|(OALL(i,41)==53)|(OALL(i,41)==54)|(OALL(i,41)==55))&(OALL(i,27)<DT)); %cancel any old day/opening/fill or kill/immediate orders 
                               %Cancellations
                            elseif...%if 1 record in 1 part and cancelled before date time in toto
                               (OALL(i,22)==1)&&(((OALL(i,40)==DT)&&(OALL(i,3)<=TM))||((OALL(i,40)<DT)))&&(OALL(i,16)>0)&&(OALL(i,30)>48)  
                                CIND(i,1)=1;
                            elseif ...%if multiple records in group and cancelled before date time
                               (OALL(i,22)>1)&&(((OALL(i,40)==DT)&&(OALL(i,3)<=TM))||((OALL(i,40)<DT)))&&(OALL(i,16)>0)  
                                CIND(i,1)=1;    
                            end

                            elseif...%Prebook records.  Pre-order date and has a tdate < event date or cancel or executed on 
                                (OALL(i,27)<901101) %Pre-901101
                                if...%%Matching order found after the date but no matching prior order
                                (isempty(find(((OALL(:,38)==OALL(i,38))&(OALL(:,30)==48)&(OALL(:,23)==OALL(i,23))&(OALL(:,31)==OALL(i,31))&(OALL(:,27)==OALL(i,27))&(OALL(:,26)==OALL(i,26))&(OALL(:,29)==OALL(i,29))&(OALL(:,32)==OALL(i,32))&(OALL(:,41)==OALL(i,41))&(((OALL(:,40)==DT)&((OALL(:,36)>TM)|(OALL(:,3)>TM)))|(OALL(:,40)>DT))), 1))==0)&...
                                (isempty(find(((OALL(:,38)==OALL(i,38))&(OALL(:,30)==48)&(OALL(:,23)==OALL(i,23))&(OALL(:,31)==OALL(i,31))&(OALL(:,27)==OALL(i,27))&(OALL(:,26)==OALL(i,26))&(OALL(:,29)==OALL(i,29))&(OALL(:,32)==OALL(i,32))&(OALL(:,41)==OALL(i,41))&(((OALL(:,40)<OALL(i,40))))), 1)))
                                PIND(i,1)=1;
                                end
                        end
                    end
                    %Duplicate orders index
                    DIND=((OALL(:,30)==48)&((OALL(:,35)==0)|(OALL(:,16)==0)));

                    %Separated partial order in 1 record.  I want to ignore these
                    SPIND=((OALL(:,35)<OALL(:,31))&(OALL(:,16)==0)&(OALL(:,22)==1)&(OALL(:,35)>0));

                    clear LB LB1 LB2 X Y
                    if BS==1
                    LB1=OALL(find((EIND==0)&(CIND==0)&(DIND==0)&(SPIND==0)&(OALL(:,27)>=901101)&(OALL(:,30)>48)&(OALL(:,23)<=QT(BS))),:);
                    LB2=OALL(find((PIND==1)&(OALL(:,23)<=QT(BS))),:);
                    LB=sortrows([LB1(:,[31 23]); LB2(:,[31 23])],[-2]);
                    %Get prices
                    clear P1 P2 Q2
                    P1=[LB(:,2)];
                    P2=P1(P1~=[1;P1(1:end-1)]);
                    Q2=NaN*ones(length(P2),1);
                    for p=1:length(P2)
                        Q2(p,1)=nansum(LB(find(LB(:,2)==P2(p)),1));
                    end
                    X=[ones(length(P2),1) log(P2) ];  %constant + log(P)
                    Y=log(cumsum(Q2)); %log(Q)
                    else
                    LB1=OALL(find((EIND==0)&(CIND==0)&(DIND==0)&(SPIND==0)&(OALL(:,27)>=901101)&(OALL(:,30)>48)&(OALL(:,23)>=QT(BS))),:);
                    LB2=OALL(find((PIND==1)&(OALL(:,23)>=QT(BS))),:);
                    LB=sortrows([LB1(:,[31 23]); LB2(:,[31 23])],[2]);
                    %Get prices
                    clear P1 P2 Q2
                    P1=[LB(:,2)];
                    P2=P1(P1~=[1;P1(1:end-1)]);
                    Q2=NaN*ones(length(P2),1);
                    for p=1:length(P2)
                        Q2(p,1)=nansum(LB(find(LB(:,2)==P2(p)),1));
                    end
                    X=[ones(length(P2),1) log(P2)];  %constant + log(P)
                    Y=log(cumsum(Q2)); %log(Q)
                    end
                    %Elasticity regression
                    clear OLSRES
                    %B=inv(X'*X)*X'*Y;
                    if length(X)>2; %require 3 observations
                    OLSRES=ols(Y,X);
                    CON(s,d,t,b)=OLSRES.beta(1);
                    E(s,d,t,b)=OLSRES.beta(2);
                    TSTATCON(s,d,t,b)=OLSRES.tstat(1);
                    TSTATE(s,d,t,b)=OLSRES.tstat(2);
                    R2(s,d,t,b)=OLSRES.rsqr;
                    N(s,d,t,b)=OLSRES.nobs;
                    BIDASK(s,d,t,b)=QT(b);
                    DEPTH(s,d,t,b)=QTSIZ(b);
                    end
                end
            end
        end
    end
end
save([TORQELASTICITYEST.mat'],'E','CON','TSTATCON','TSTATE','R2','N','BIDASK','DEPTH')
